subroutine mt_cohe1 ( npts,dt,xi,xj,tbp,kspec,nf,p,                   &
                     freq,spec,cohe,phase,speci,specj,conf,               &
                     cohe_ci, phase_ci,iadapt )  
!
!  Construct the coherence spectrum from the yk's and the 
!  weights of the usual multitaper spectrum estimation. 
!  Note this code uses the real(4) multitaper code. 
!
!  INPUT
!
!	npts		integer number of points in time series
! 	dt		real, sampling rate of time series
!	xi(npts)	real, data for first series
!       xj(npts) 	real, data for second series
!	tbp		the time-bandwidth product
!	kspec		integer, number of tapers to use
!	nf		integer, number of freq points in spectrum
!	p		confidence for null hypothesis test
!
!
!  OPTIONAL INPUT
!
!	iadapt		integer 0 - adaptive, 1 - constant weights
!			default adapt = 1
!
!  OPTIONAL OUTPUTS
!
!	freq(nf)	real vector with frequency bins
!       cohe(nf)	real, coherence of the two series (0 - 1)
!       phase(nf)	the phase at each frequency
!	speci(nf)	real vector with spectrum of first series
!	specj(nf)	real vector with spectrum of second
!	conf(nf)	p confidence value for each freq.
!	cohe_ci(nf,2)	95% bounds on coherence (not larger than 1)
!	phase_ci(nf,2)  95% bounds on phase estimates
!
!	If confidence intervals are requested, then both phase and
!	cohe variables need to be requested as well. 
!
!  Modified
!
!	German Prieto
!	September 2005
!
!	*******************************************************************
!
!	German Prieto
!	October 2007
!  
!	Re-wrote the subroutine to allow optional output argu-
!	ments. 
!
!	******************************************************************* 
!
!	German Prieto
!	February 2008
!  
!	Added optional adaptive or constant weighting, similar to 
!       mt_transfer.
!
!	******************************************************************* 
!

!  calls
!	mtspec
!	


!********************************************************************

   use spectra

   implicit none

!  Inputs

   integer, intent(in) :: npts, kspec, nf

   real(4), intent(in) :: dt, tbp
   real(4), intent(in out) :: p

   real(4), dimension(npts), intent(in) :: xi, xj

!  Optional Input

   integer, optional                        :: iadapt

!  Optional outputs

   real(4), dimension(nf), optional :: freq, speci, specj

   real(4), dimension(nf), optional :: cohe, phase, conf
 
   real(4), dimension(nf,2), optional :: cohe_ci, phase_ci

!  spectra and frequency

   real(4), dimension(nf) :: f, si, sj

   real(4), dimension(nf,kspec)      :: wt_i, wt_j 

   complex(4), dimension(npts,kspec) :: yk_i, yk_j

   real(4), dimension(2)             :: xmean, xvar

   real(4), dimension(nf)            :: wt_scale


!  Coherence freq matrices

   complex(4), dimension(nf,kspec)     :: dyk_i, dyk_j 
   complex(4), dimension(nf), optional :: spec

!  Confidence terms

   real(4), dimension(kspec)         :: wt_dofs
   real(4)                           :: se

!  Jackknife

   complex(4), dimension(nf,kspec) :: cross
   real(4), dimension(nf,kspec)    :: spec1, spec2

   complex(4), dimension(nf)       :: cross_all
   real(4), dimension(nf)          :: spec1_all, spec2_all
   real(4), dimension(nf)          :: cohe_all, phase_all

   real(4), dimension(nf,kspec)    :: cohejk, Q 
   real(4), dimension(nf,kspec)    :: qvarjk
   real(4), dimension(nf)          :: qmean, qvar, phvar
   real(4), dimension(nf)          :: qtop, qbot, qval

   complex(4), dimension(nf,kspec) :: phasejk
   complex(4), dimension(nf)       :: phmean

   real(4) :: atanh2
   
!  Others

   integer :: i, j, iad

   real(4), parameter :: pi = 3.14159265358979

!********************************************************************

!   write(6,'(a)') 'Calling multitaper coherence code'
!   write(6,'(a,2i6)') 'Data points and frequency points', npts, nf

   if (present(cohe_ci) .and. .not. present(cohe)) then
      write(6,'(a)') 'If cohe_ci is requested, cohe is needed, stopped'
      stop
   elseif (present(phase_ci) .and. .not. present(phase)) then
      write(6,'(a)') 'If phase_ci is requested, phase is needed, stopped'
      stop
   endif

   if (p>=1. .or. p<=0.) then
      p = 0.95
   endif
!   write(6,'(a,f10.7)') 'Confidence for null hypothesis testing ', p


   xmean(1) = sum(xi)/real(npts)
   xvar(1) = (sum((xi - xmean(1))**2))/real(npts-1)
   xmean(2) = sum(xj)/real(npts)
   xvar(2) = (sum((xj - xmean(2))**2))/real(npts-1)

! Adaptive or constant weighting?

   if (present(iadapt)) then 
      if (iadapt == 0) then
         iad = 0
      else
         iad = 1
      endif
   else
      iad = 1
   endif

!
!  Get the spectrum estimate
!

!   if (all(xi==xj)) then    
!      write(6,'(a)') 'Auto-coherence not applicable'
!      stop
!   else
      call mtspec_r( npts,dt,xi,tbp,kspec,nf,f,          &
                    si,yk=yk_i,wt=wt_i,adapt=iad)

      call mtspec_r( npts,dt,xj,tbp,kspec,nf,f,          &
    		    sj,yk=yk_j,wt=wt_j,adapt=iad)
!   endif

!
!  Create the spectra (cannot use spec output, normalized different)
!

   wt_i = min(wt_i,wt_j)
   wt_j = min(wt_i,wt_j)

   wt_scale = sum(wt_i**2, dim=2)  ! Scale weights to keep power 
   do i = 1,kspec
      wt_i(:,i) = wt_i(:,i)/sqrt(wt_scale)
      wt_j(:,i) = wt_j(:,i)/sqrt(wt_scale)
   enddo

   do i = 1,nf
      do j = 1,kspec
         dyk_i(i,j) = wt_i(i,j) * yk_i(i,j)
         dyk_j(i,j) = wt_j(i,j) * yk_j(i,j)
      enddo
   enddo

   si = sum(abs(dyk_i)**2, dim=2) 
   sj = sum(abs(dyk_j)**2, dim=2) 

      freq = f
   do i=1,nf

   ! coherence and phase
   
      spec(i) = sum ( dyk_i(i,:) * conjg(dyk_j(i,:)) )  


   enddo
end subroutine mt_cohe1

